Preparations

Load the necessary libraries

library(car)       #for regression diagnostics
library(broom)     #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot)    #for outputs
library(knitr)     #for kable
library(effects)   #for partial effects plots
library(ggeffects) #for partial effects plots
library(emmeans)   #for estimating marginal means
library(MASS)      #for glm.nb
library(MuMIn)     #for AICc
library(tidyverse) #for data wrangling
library(nlme)
library(lme4)
library(glmmTMB)
library(broom.mixed)
library(glmmTMB)   #for glmmTMB
library(DHARMa)   #for residuals and diagnostics
library(performance) #for diagnostic plots
library(see)         #for diagnostic plots

Scenario

To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.

A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.

Barramundi

Sampling design

Format of norin.csv data files

FISHID MASS TRIAL SMR_contr CHANGE
1 35.69 LowSalinity 5.85 -31.92
2 33.84 LowSalinity 6.53 2.52
3 37.78 LowSalinity 5.66 -6.28
.. .. .. .. ..
1 36.80 HighTemperature 5.85 18.32
2 34.98 HighTemperature 6.53 19.06
3 38.38 HighTemperature 5.66 19.03
.. .. .. .. ..
1 45.06 Hypoxia 5.85 -18.61
2 43.51 Hypoxia 6.53 -5.37
3 45.11 Hypoxia 5.66 -13.95
FISHID Categorical listing of the individual fish that are repeatedly sampled
MASS Mass (g) of barramundi. Covariate in analysis
TRIAL Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen.
SMR_contr Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia)
CHANGE Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'.

Read in the data

norin = read_csv('../public/data/norin.csv', trim_ws=TRUE)
glimpse(norin)
## Rows: 180
## Columns: 5
## $ FISHID    <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
## $ MASS      <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.…
## $ TRIAL     <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity",…
## $ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589,…
## $ CHANGE    <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.…

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of temperature and (centered) mean fish size on SDA peak. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual fish.

Fit the model

For each of the following modelling alternatives, we will:

  • establish whether MASS is a useful covariate (based on AICc). This step involves evaluating alternative fixed effects structures. When doing so, we must use Maximum Likelihood (ML) rather than Residual Maximum Likelihood (REML).
  • establish whether a random intercept/slope model is useful (based on AICc). This step involves evaluating alternative random effects structures. When doing so, we should use Restricted Maximum Likelihood (REML)

lme (nlme)

fixed structures

##Compare models that estimate partial slope for MASS vs an offset for MASS
##must use ML to compare models that vary in fixed effects
norin.lme1 <- lme(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + offset(MASS), random=~1|FISHID, data=norin, method='ML')
## Unfortunately,  update() does not remove offset().
## We will just have to write the other model out in full as well.
norin.lme2 <- lme(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + MASS,  random=~1|FISHID, data=norin, method='ML')
## Now without MASS altogether
norin.lme3 <- lme(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE),  random=~1|FISHID, data=norin, method='ML')

## Compare these models via AICc
AICc(norin.lme1,norin.lme2, norin.lme3)
## Alternatively,  we can use sequential Likelihood Ratio Tests (LRT)
anova(norin.lme1, norin.lme2, norin.lme3)

Conclusions:

  • on the basis of AICc, the model without MASS is considered the most parsimonious (lowest AICc)
  • we will proceed with the fixed structure of model 3

random structures

norin.lme3a = update(norin.lme3, method='REML')
norin.lme3b = update(norin.lme3a, random=~TRIAL|FISHID)
AICc(norin.lme3a,  norin.lme3b)
anova(norin.lme3a, norin.lme3b)

Conclusions:

  • the AICc of the random intercept/slope model is lowest
  • therefore we will proceed with a random intercept/slope model

lmer (lme4)

fixed structure

##Compare models that estimate partial slope for MASS vs an offset for MASS
##must use ML to compare models that vary in fixed effects
norin.lmer1 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + offset(MASS) + (1|FISHID), data=norin, REML=FALSE)
## Unfortunately,  update() does not remove offset().
## We will just have to write the other model out in full as well.
norin.lmer2 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + MASS + (1|FISHID), data=norin, REML=FALSE)
## Now without MASS altogether
norin.lmer3 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE)+ (1|FISHID), data=norin, REML=FALSE)

## Compare these models via AICc
AICc(norin.lmer1,norin.lmer2, norin.lmer3)
## Alternatively,  we can use sequential Likelihood Ratio Tests (LRT)
anova(norin.lmer1, norin.lmer2, norin.lmer3)

Conclusions:

  • on the basis of AICc, the model without MASS is considered the most parsimonious (lowest AICc)
  • we will proceed with the fixed structure of model 3

random structure

norin.lmer3a = update(norin.lmer3, REML=TRUE)
## unfortunately,  the following random intercept/slope model cannot be fit
##norin.lmer3b = update(norin.lmer3a, ~ . -(1|FISHID) +(TRIAL|FISHID))
##AICc(norin.lmer3a,  norin.lmer3b)
##anova(norin.lmer3a, norin.lmer3b)

Conclusions:

  • unfortunately, we are unable to fit a model with random intercept/slope
  • therefore we will proceed with a random intercept model

glmmTMB

fixed structure

##Compare models that estimate partial slope for MASS vs an offset for MASS
norin.glmmTMB1 = glmmTMB(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + offset(MASS) + (1|FISHID), data=norin, REML=FALSE)
## Unfortunately,  update() does not remove offset().
## We will just have to write the other model out in full as well.
norin.glmmTMB2 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + MASS + (1|FISHID), data=norin, REML=FALSE)
## Now without MASS altogether
norin.glmmTMB3 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + (1|FISHID), data=norin, REML=FALSE)

## Compare these models via AICc
AICc(norin.glmmTMB1,norin.glmmTMB2, norin.glmmTMB3)

Conclusions:

  • on the basis of AICc, the model without MASS is considered the most parsimonious (lowest AICc)
  • we will proceed with the fixed structure of model 3

random structure

norin.glmmTMB3a <- update(norin.glmmTMB3, REML=TRUE)
## unfortunately,  the following random intercept/slope model cannot be fit
## norin.glmmTMB3b <- update(norin.glmmTMB3a, ~ . - (1|FISHID)+ (TRIAL|FISHID))
##AICc(norin.glmmTMB1a,  norin.glmmTMB1b)

Conclusions:

  • unfortunately, we are unable to fit a model with random intercept/slope
  • therefore we will proceed with a random intercept model

Model validation

Partial plots model

Model investigation / hypothesis testing

Predictions / further analyses

Summary figures

References

Norin, Tommy, Hans Malte, and Timothy D. Clark. 2015. “Differential Plasticity of Metabolic Rate Phenotypes in a Tropical Fish Facing Environmental Change.” Functional Ecology 30 (3): 369–78. https://doi.org/10.1111/1365-2435.12503.